function [ fout, tout ] = reconstruction( fin, fsin, fsout )

Tout = 1 / fsout;
Tin = 1 / fsin;

if (fsin == fsout)
    fout = fin;
    tout = Tout* linspace(0, size(fout, 2), size(fout, 2)); 
    return;
end

scale = fsout / fsin;

if (scale < 1)
    % for compatibility with downsampling and decimate
    fout = zeros(1, ceil(size(fin, 2)* scale));
else
    fout = zeros(1, fix(size(fin, 2) * scale));
end

tout = linspace(0, (1/fsin) * (size(fin, 2) - 1), size(fout, 2));            

% when downsampling cut frequencies below fsout/2 from fin 
% function.
if (fsin > fsout)
    k = sinc_gen(fsin, fsout/2);
    fin = conv(fin, k, 'same');
end

W = fix(160 / (fsout*2 / fsin)); % experimental
W = min(W, size(fin, 2));

tin = (0:1:size(fin, 2)-1) * Tin;

if (W > size(fin, 2))
    W = size(fin, 2);
end

for m=1:size(fout, 2)
    [~, min_index] = min(abs(tout(m) - tin));
    min_index = fix(min_index - W/2);
    if (min_index < 1)
        min_index = 1;
    end
    max_index = min_index + W - 1;
    if (max_index > size(fin, 2))
        max_index = size(fin, 2);
    end
    sc = sinc((tout(m) - (min_index-1:1:max_index-1)*Tin) / Tin);
    fout(m) = sum(fin(min_index:1:max_index) .* sc);
end

% for compatibility with upsampling
if (scale > 1)
    fout = fout(1:1:end-scale+1);
    tout = tout(1:1:end-scale+1);
end

